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Abstract — The state of art of time domain integral equation 
(TDIE) solvers has grown by leaps and bounds over the past 
decade. During this time, advances have been made in (i) the 
development of accelerators that can be retrofitted with these 
solvers and (ii) understanding the stability properties of the 
electric field integral equation. As is well known, time domain 
electric field integral equation solvers have been notoriously 
difficult to stabilize. Research into methods for understanding and 
prescribing remedies have been on the uptick. The most recent of 
these efforts are (i) Lubich quadrature and (ii) exact integration. 
In this paper, we re-examine the solution to this equation using 
(i) the undifferentiated form of the TD-EFIE and (ii) a separable 
approximation to the spatio-temporal convolution. The proposed 
scheme can be constructed such that the spatial integrand over 
the source and observer domains is smooth and integrable. As 
several numerical results will demonstrate, the proposed scheme 
yields stable results for long simulation times and a variety of 
targets, both of which have proven extremely challenging in the 
past. 

Index Terms — Time Domain Analysis, Integral Equations, 
Marching on in Time, Separable Expansions, Stability. 



I. Introduction 

Since the development of TDIE based methods in the 
1960's [1], computational complexity and late time instability 
have been the two stumbling blocks that have prevented their 
widespread adoption as electromagnetic analysis tools. Over 
the past decade, the former has been largely addressed and fast 
evaluators with provable error estimates that can be retrofitted 
with TDIE solvers exist [2,3]. In fact, TDIE solvers aug- 
mented with these fast solvers have seen widespread applica- 
tion to a range of electromagnetics and acoustics applications 
[4, 5]. While the issue of computational complexity is a solved 
problem, that of instability still lingers. Over the years there 
have been several attempts at stabilizing these equations; 
these include filtering techniques [6], implicit time stepping 
[7], smooth basis functions [8], space-time Galerkin methods 
[9], and Lubich quadrature [10]. Methods that are based on 
transform techniques have been developed as well [7], and they 
avoid problems with instability by avoiding time marching 
altogether. As they avoid marching on in time (MOT), they 
will not be the subject of the ensuing discussion. 
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Despite the plethora of interest, instability of TDIE solvers 
is far from being a solved problem. Of all the methods listed 
above, only the space-time Galerkin methods are provably 
stable [11]. Here, stability is proven by passage through the 
Fourier-Laplace domain. Indeed, most of the others, with the 
possible exception of Lubich quadrature, simply delay the 
onset of instability. While the stability properties of space-time 
Galerkin based schemes have been proven, such methods are 
notoriously difficult to implement. These methods rely on two 
aspects to obtain stability; (i) construction of the TDIE based 
on a variational formulation and (ii) exact evaluation of all 
integrals involved. The latter involves the exact evaluation of 
five-dimensional integrals [9]. However, as most of literature 
on this topic is published only as Ph.D theses, only sparse de- 
tails have been available to the rest of the scientific community. 
Recent papers have sought to address some of these concerns 
[11], but are largely restricted to problems in acoustics. In 
the electromagnetics community, attempts to develop such a 
scheme started with the development of space-time collocation 
together with exact evaluation of fields for both scattering from 
perfectly conducting and composite [12] objects. It was shown 
that this method was indeed stable for a whole class of targets 
that had, in our experience, been impossible to stabilize with 
any of the other existing techniques. More recently, a space- 
time Galerkin formulation for electromagnetics was introduced 
[13], wherein stability was shown for extremely long solution 
times, small time steps and higher order temporal interpolants. 

The key in implementing this method is determining the 
topology of the domain of integration. This is best illustrated 
by considering the scheme within a collocation framework as 
done in [12]. Here, to test the field received by a triangular 
patch due to a point source whose time signature is piecewise 
continuous, one has to find the corresponding domains in 
the triangle where the integrands are piecewise continuous. 
This is tantamount to finding intersections of the triangle with 
concentric time spheres of radii icA t that are centered at the 
source point, where c is the speed of light, A t is time step 
size, and i is an integer. It is readily apparent that while this 
method is very effective, it is impossible to use in a higher 
order framework and extremely difficult to integrate with fast 
methods. Another approach that has seen recent attention is the 
application of Lubich quadrature [10] to TDIEs. This method 
relies on a series of transforms, to the Laplace domain, the 
z-domain, and back into the time domain. As it is based on 
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a series of domain transformations, it is possible to analyze 
propagation through and scattering from dispersive and lossy 
materials. On the flip side, as it is based on entire-domain 
transforms, the method converts what should be an 0(N t N%) 
solver given the compact nature of the retarded potential to a 
system whose complexity scales as 0{N^N 2 s ). Here, N t and N s 
are the number of temporal and spatial degrees of freedom, 
respectively. To overcome this bottleneck, one would then have 
to take recourse to an FFT based method [14]. 

The goal of this paper is to present a method that yields 
stable results for the time domain electric field integral equa- 
tion (TDEFIE) while obviating some of the aforementioned 
drawbacks. It is based on the undifferentiated form of the 
TDEFIE (and to a large extent inspired by [12, 13]), and relies 
on deriving an alternate representation of convolution between 
the retarded potential and the space-time basis function. Key 
attributes of this method are the following: (i) it does not 
transform the problem to the Laplace domain, (ii) it has finite 
support in time, and (iii) it is entirely numerical, and as a 
result, it can be extended to higher order discretizations (in 
space and time) and integrated with existing fast methods. The 
principal contributions of this paper are threefold: 

1) We will present a methodology based on a separable 
approximation to the space-time convolution of a source 
with the retarded potential. We will elucidate its imple- 
mentation within an MOT framework. 

2) We will frame the MOT system as an eigenvalue prob- 
lem (akin to that done in the [15] but for the system of 
equations used here). This will be used to provide insight 
into the stability of the system for a given discretization 
and time step size. 

3) Finally, to demonstrate the effectiveness of the method, 
we will present scattering results from a number of 
benchmark targets, with RCS comparisons to a validated 
frequency domain solver or analytical results where 
possible. A number of the presented targets have defied 
all stabilization schemes apart from [10, 12, 13]. 

In this paper, we shall not present results of extension of this 
method to either higher order surfaces or integration with fast 
solvers. Likewise, we will rely on demonstration of stability 
via a set of challenging targets as opposed to mathematical 
proofs for late time stability. 

The rest of this paper is organized as follows: Section 
II details the approach used in this paper, contrasts it with 
methods used elsewhere. Implementation details of this ap- 
proach are presented in Section II-B. Details of the eigenvalue 
analysis implementation is presented in Section III. Section 
IV presents a number of results that demonstrate late time 
stability. Lastly, Section V summarizes the contribution of this 
paper and outlines directions of future research. 

II. Formulation 

Consider a perfectly electric conducting body occupying 
a domain D- c R 3 in free space (jjq,so). Let D denote 
the bounding surface of £>_ and let the outward pointing 
unit vector normal at any point r on Q be denoted by h(r). 
Assume that a field {E*(r, t), H'(r, t)} incident on this body is 



bandlimited to frequency f max and is zero for t < 0. This field 
induces currents, J(r, t), on the surface D, that radiate scattered 
fields {E*(r, t), EP(r, t)}. Thus, the total electric and magnetic 
fields in all of D + = R 3 /Z)_, denoted by E(r,f) and H(r,f), 
can be decomposed into the incident and scattered fields as 



E(r,0 = E<(r,0 + E*(r,0 , 
H(r, t) = H (r, t) + H*(r, t) . 



(la) 



The currents induced can be solved for using either the electric, 
magnetic or combined field integral equation. The TDEFIE has 
historically been the most challenging to stabilize. Therefore, 
the rest of this paper will focus on discretizing this equation. 
The TDEFIE may be written as 



{E ,- (r, t), H (r, t)} 




{E 5 (r,0, H 5 (r,f)} 



Fig. 1: General description of a scattering problem 



h x hxE\r, t) = -h x h x E s o {J(r, t)} Vr e Q , 
E* o {J(r, 0} = ~d t A o { J(r, 0} - VO o { J(r, 0} , 



If f T V • 

Oo{J(r,f)}= dr' dt' 

47T£ Jn J-oo 



(lb) 



J(r',0 



R 



where R = |R| = |r-r'|, r = t-R/c and V denotes a divergence 
with respect to r' . The solution to this integral equation for the 
unknown currents J(r, t) is typically effected by representing 
the current in terms of space-time basis functions as 



N t 



J(r,0 = X S " (r) Z 7 ^ (0 ' 



(2) 



where the temporal coefficients J n j are to be determined for 
N t temporal and N s spatial basis functions. In the above 
expressions, the spatial basis functions S n (r) are chosen to be 
the Rao-Wilton-Glisson basis functions [16] defined for each 
edge on the tesselation that describes Q such that 

— (r - r + ) reP + 

— — (r - r") reP" 



S„(r) = 



(3) 



where P„ are the two triangles (of areas A*, respectively) 
that are associated with each edge n, l n is the length of the 
edge, and are free vertices associated with P*. The spatial 
basis function vanishes outside the domain £l n = P+UP~. The 
temporal basis functions used, T t {t) = T(t - iA t ), are typically 
pth order shifted Lagrange polynomials given by [12] 



T{t) =fk(t)g P -k(t)P a ,p(t) , 

a = (k- V)A t \P = kA t for k = 0, • 



(4) 
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where P a ^(t) is a rectangular pulse function in the domain 

1 k = 



fk(f) 



7=1 V 



t ~ jA t 



k±0 ' 



and 



P -k 

gp-k = I J 

7=1 



(5) 



(6) 



where A t = x/(20fmax) is the time step size and x i s an 
oversampling factor. Traditional MOT schemes are derived 
by substituting (2) in (lb), using Galerkin testing in space 
and point testing in time. The resulting equations may be 
succinctly written as 

7-1 7-1 

zo i j =fj-Yj ZiJ H - Z tiC ^ ' (7a) 



i=\ 



where 



1 j ~ [Jlj'hj,' - >Jn s ,j] , 

rJ At 



(7b) 



:j = cj-i + v it dt'uf) , (7c) 

r nJ = (S n (r),axaxE f (r,0>| , (7d) 

\ / \t=jA t 



and 



Z** = (S w (r),iix h x EJ o {S w (r)r^.(0})| ' 
Znw,,- = (S n (r), nxnxEJo {s m (r)T>_ f (0})| ' 



(7e) 



E 5 = E* + E* , 
where E* 2 are defined as 

E^ o {s w r f } = - d t A m o {s w r f } - voi o {s w r f } , 
e* o {s w r f } = - vo 2 o {s w r f } , 

<*i o {Sm r,} = J- f ^ v/ ' s ; (r/) P at^ ) , (7f) 

, /C1 1 f , V'-S m (r') 

$2°{S m }=- dr , 

47i£ Jq R 

where A: = |r/AJ. In the above equations, (•) denotes a 
standard inner product. As is evident from (7a), the solution 
proceeds sequentially over each time step. The crux of the so- 
lution of late time instability has been the accurate evaluation 
of (7e) (together with the use of the undifferentiated TDEFIE). 
It is evident from the nature of the temporal basis functions 
that the integrand in (7e) is piecewise-continuous, and the 
use of Gauss quadrature to evaluate these integrals before 
first identifying domains of continuity will be inaccurate. It 
should be noted that the situation is far worse when one uses 
the derivative form of the TDEFIE as the vector potential 
term includes the second derivative of the temporal basis. In 
what follows, we will seek to develop an alternate method for 
evaluating this integral. 



A. Approximation of the temporal convolution 

The approach espoused in this paper is to approximate the 
convolution of the space time basis function with the retarded 
potential. Note, that the TDEFIE in (lb) requires both the 
convolution with the temporal derivative and the temporal 
integral of this space-time basis, and expressions for effecting 
these will be provided as we proceed. For simplicity, consider 
a point source and an observation triangle as shown in Fig. 2. 
Consider a field due to the point source given by 



<A(r, = 



S(t-R/c) 
4nR 



** T t (t) , 



(8) 



where * r denotes a convolution with respect to time. Given 
the piecewise nature of T t (t), it is evident that the lines of 
discontinuity correspond to intersections of time spheres with 
the observation triangle. This scheme has been implemented 
in acoustics [11, 17, 18] and in elecromagnetics [12] to great 
success. 

The approach proposed herein takes a slightly different path. 
Define the radii of the smallest and largest time spheres that 
enclose the triangle as (acA t + f ) and (j3cA t + £), respectively, 
where f is the largest multiple of cAt between the source 
point and observer triangle. Within this region, the convolution 
with the retarded potential can be expressed as a separable 
expansion in space and time. Using this expansion, the field 
due to a point source located at r' with temporal dependence 
Ti(f) may be approximated as 

AnR 



AnR 



1=0 



i4 _ 9*'I> P/(fw)f < (0 ' 

1=0 



(9) 



where i(R) = k\ {R - £)/c + k 2 , a\-k\ 



2/+ 1 



t\(t) = Ptiht + k 2 yp a f(t/A t ) ut) . 

The constants k\ and k 2 are chosen such that they map the 
domain [a,/3] — > [-A t ,A t ], Pi(-) is a Legendre polynomial of 
order / and Nh is the number of harmonics that are retained 
in the expansion. The key consequence of this expression 
is the separation of space and time within the domain R e 
[acA t + £,/3cA t + £]. As will become evident, this leads to 
spatial integrands which are smooth over the entire triangle. 
It is evident that this can be generalized to the case of a 
source/observation pair. For instance, 



A(r,f)o{S w r f } = 



*(,-M) 



47T 

A- 1) 



* st S m (r)Ti(t) , 



-£i>K) 
*(>-¥) 



R 



S m (r')Ti(t) , 



* t S m (r')Ti(t) . 



(10) 
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Using (9) in the last term, it follows that 

Nh 

<S„(r), A(r, t) o {S m m « 5 (f - |) g g a^ffr) , (1 la) 
where 

Here f , a and /3 are chosen such that < (R- £) /c < f3A t 
for all r' e Q m and r e From the above expression 



where 



Source 







R-i \ 


X x\ 

/ 1 


cAt 





Observation 



Fig. 2: Arcs of intersection 



it is evident that the spatial and temporal convolutions are 
completely independent of each other. It is also observed 
that the spatial convolution is smooth over the observation 
domain with a removable singularity. Given that one a priori 
knows Nh, integration rules can be designed to evaluate 
to high accuracy. The piecewise continuity or discontinuity 
of the temporal basis functions is present only in temporal 
convolutions and is handled relatively easily. 



B. Implementation Details 

The above exposition details evaluation of the tested vector 
potential. Extensions to evaluate the specific components of 
(7a) are prescribed next. The matrix elements in (7a) comprise 
the temporal derivative of the vector potential and the scalar 
potential. Specifically, using (11a), it follows that 



(12) 



t=jA t 



Note, the temporal basis function is compact. This convolution 
(at integer multiples of A t ) can be evaluated analytically. The 
evaluation of the contribution due to the scalar potential is a 
little more involved. To this end, contribution of the scalar 
potential can be written as 

<S B (r), VO(r, t) o {S m Ti})\ t=jAt =s[t-^* t 



> a &Pi(ht + k 2 yP a fi(t/At)* t (13a) 



J -ex 



dt'm') 



-f 



drW • S n (r) 



1 



dr 



R 



(13b) 



Here the vector derivative in r has been moved onto the testing 
function. The evaluation of the convolution of the Legendre 
polynomial with the integral of the basis function may be 
written as 



PiQdt + k 2 )P afi (f/A t ) *, f dt'Ti(t') = 

*J — oo 

Xoo W-T 
dTP l (k l T + k 2 )P afi (T/A t ) df 
CO U —CO 



(14) 



Ti(t') 



It can be shown that 



r 

*J —oo 



dt'Ti(t') = 



dt'Tiit') t > yA, 

oo 

| dt'Ti(t') t < yA t 

^ \J —CO 



(15) 



where y = ft - a + i + p. As a result, the integrals in (14) 
decouple into a product of two integrals for t > yA t , and can 
be evaluated using the fact that drPi(k\T + kifP a fi{j j 'A t ) = 
doj(fi-a)At, where 6ij is the Kronecker delta function. Using 
this, the values for sufficiently large t can be moved to the 
right hand side of (7a), analogous to the procedure in (If). 
The relation (7c) becomes 

Cj= Cj-x+Ij-x dt f T hl (t f ) (16) 

J(j-2)At 

Finally, as an aside, we note a couple of other implementation 
features: 

1) The temporal convolutions, while associated with the 
source and testing function, are independent of space. 
Furthermore, these are always tested using point testing 
at integer time steps. This implies that they need to be 
evaluated only once for each value of p - a, and can be 
precomputed for a given geometry. 

2) The highest order of the polynomial is known for a 
pair of basis functions and is given by A^. As a result, 
rules for evaluating these integrals can be obtained from 
standard libraries or constructed. Furthermore, for self- 
triangles, the sinh -1 rule outlined in [19] is used. Exten- 
sion of this method to the integrand of each harmonic 
is trivial and is not elaborated here. 

III. Eigen Analysis of the MOT system 

In order to analyze the stability of the scheme proposed 
here, we perform an Eigen spectrum analysis similar to that 
in [15], but for the form of equations given here. Consider 
a simplification of the marching scheme in equation (7a) and 
modified with the separable expansion, given by 

7-1 7-1 

zoij =fj~Yj ZiI ^ - Z tiC ^ ' (17a) 

i= 1 i= 1 

where Cj are the coefficients of the charge at time step j and 
Zi has been modified from that in (If) using the relation in 
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(15). From the equation of continuity we have that the charge, 
p(r, ti), is given by 



p(r, ti) 



= f dt'V-J(r,f) , 

%J — oo 

Xt~At r>ti 
dt'V -J(r,r')+ dt'V- J(r,r') , 
oo Jti-At 



(17b) 



= p(r, - Af) + 



J^V- J(r,0 , 



which using the basis function expansions can be expressed as 
a relationship between the temporal coefficients of the charges 
and currents as defined in (16). To start our analysis, we first 
note that the summations on the right hand side of equation 
(17a) can be restricted to run from / = 1 to / = j - P, where 
P = max(l, N max ), where N max , in turn, is the temporal extent 
of the scatterer. In other words, (N max - l)cAt < diam(Q) < 
(N max )cAt, where diam(Q) is the spatial extent of the scatterer 
along the incident field direction. Then equations (17a) and 
(16) can be written in matrix form as 



" An 







Bu 


B\2 




An 


I 


h = Zi- 


B21 


§22 


Ij-i , (18) 



where An, A 2 \, B n , B21, B u , B 22 , I and !F/ are defined 
in the appendix. We denote the two matrices on either side of 
equation (18) as A and B and assume that the temporal vector 
Ij contains both current and charge coefficients. Assuming 
that the matrix A is invertible, and denoting C = A~ l B, the 
current vector at any given time step can be written as 

I_x = A-^-CTo, 

12 = A- x n-c{A- x n-ci^ , 

j 3 = A~ l r 3 - cA~ l r 2 + c 2 A~ l ri - c 3 i , (19) 



7-1 

Now, let the eigenvalue decomposition of C be given by [20] 
C = ^V^, (20a) 



where t represents a conjugate transpose. Then, 



(20b) 



can be used to provide a natural bound on the matrix vector 
product as 

ll^dM^Ev^l (20c) 



for any vector ( F, where cr is the largest eigenvalue of C. 
Using the bound in (20c), Ij is bounded by 

7-1 



k=0 q 



(21) 



Assuming that after some time jA t the input signal vanishes, 
i.e., Tj - 0, then the first summation on the right hand side of 
equation (21) can be restricted to k e p\,p\ + 1, . . . ,p2 where 
the temporal extent of the input signal is from p\ to p 2 . Then 
(21) leads to 



:i + ^ 



(22) 



where P = p 2 - p\ + 1 and Q = v q ^v q max^ {a -1 ^} J. At 
this point, it is apparent that depending on the value of cr , 
any numerical error in the initial current Jo will grow without 
bound if cr > 1.0, will decay to a value strictly bounded 
above by PCi is <t < 1.0, and will have a constant DC value 
bounded above by PCi if cr = 1. 

IV. Results 

In this section, we present a number of results of scattering 
by perfect electrically conducting objects that are topologically 
very different from each other, some of which have, to the 
best of our knowledge, defied all stabilization attempts in 
the past with the exceptions of Lubich quadrature and exact 
integration or its variation [10, 12, 13]. Three pertinent features 
of the results presented are as follows: (i) the excitations 
are broadband, (ii) the solution time is very long (enough to 
permit multiple transits across the object), and (iii) some of the 
scatterers analyzed have features that make analysis difficult. 
The singular purpose of the objects chosen for analysis is to 
demonstrate late time stability for multiple transits on very 
challenging targets that are excited by a broad band pulse. In 
all cases, the incident field is a plane wave of the form 

E f (r, t) = u cos(27r/o0e~ ( '~ r '* /c ~^ )2/2£r2 , (23) 

where u denotes the polarization vector, / the center fre- 
quency, k = sin(^) cos(0')i; + sin(^) sin(0*)j + cos(6 l )z the 
direction of propagation, and l and (p l the polar and azimuthal 
angles of incidence, respectively. The values cr and t p are 
calculated as cr = 3/(2nB) and t p = 6<r, where B denotes 
the bandwidth of in Hz. The incident power is calculated 
to be approximately 160 dB below the peak at f max = fo + B 
and f m i n = fo - B. For all results obtained in this paper, first 
order Lagrange polynomials are used as basis functions. In 
each of the examples, we present the current at a point on the 
geometry as a function of time. In addition, for some of the 
cases analyzed, we present an eigen spectra that shows that all 
values lie within the unit circle, with the exception of few that 
approach 1.0. Finally, for some of the benchmark targets, we 
extract frequency domain radar scattering cross-section (RCS) 
data at three different frequencies, and this is then compared 
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with similar data obtained either from an analytical code or a 
validated frequency domain integral equation solver. 

A. Sphere 

Our first test case is a sphere of radius 1 m discretized with 
576 unknowns. It is excited with an incident wave propagating 
along k = z, polarized along u = x, with f min = 10 kHz, 
and f max = 182 MHz. Shown in Fig. 3 is the current density 
at (x,y,z) = (1,1.29,1.70) m for 47,000 time steps, with 
X = 4/3. This example is also run for a much larger time step 
corresponding to x - 4. The inset shows the same current 
zoomed in the early time until .45 yus. As is evident from 
this figure, the currents are stable, expectedly so for the larger 
time step, as well as the smaller time step. The two currents 
agree well with each other. Fig. 3 shows that the surface 




% s) 

Fig. 3: Current on 1 m radius sphere 

current remains bounded over the duration of the simulation. 
To show that the result is indeed stable, we performed the 
analysis outlined in section III. Fig. 5 shows that all of the 
eigenvalues lie within the unit circle and as a result, stability 
can be expected. Next, Fig. 6 compares the RCS of the sphere, 



1 




-1 -0.5 0.5 

Fig. 4: Eigenvalues of sphere MOT matrix 



observed in the x - z plane for 6 e [-180°, 180°], against that 
obtained using Mie series for three different frequencies. As is 
evident, agreement between the two sets of data is excellent. 
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Fig. 5: Convergence of RCS of sphere and plate with respect 
to* 

Fig. 5 shows the convergence with respect to oversampling 
factor*, of the RCS of the sphere at / = /o as well as a lxl 
m plate discretized with 133 unknowns at / = fo = 75 MHz. 
Convergence here is defined as Q = \\cr(Xi)-cr(Xi-i)\\/\\cr(Xi)\\ 
where, ||-|| is the norm and RCS = 10/<9gio|<x|. No more than 
20 harmonics were used in any source/observer interactions in 
these simulations. 




O Analytical ; f = 40 MHz 
x 90 MHz 
* 140 MHz 
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Fig. 6: RCS of 1 m radius sphere 



B. Two Plates 

Next, we apply this scheme to analyze scattering from an 
object described using an open surface. Here, the object chosen 
was similar to that in [13] and consists of two plates of size 
lm x lm parallel to the x - y plane that are separated by 
0.1m. Each plate is discretized using 560 spatial unknowns. 
The incident field is polarized along u = \{x + y - y/lz), 
propagates along k = ^(x + y + ^flz), and is described by 
parameters f min = 100 kHz and f max = 264 MHz. The current 
is observed at (x,y,z) = (0.95,0.95,0.1) m for 100,000 time 
steps for a time step given by x = 1- This corresponds to 
approximately 3,383 transits across the object after the incident 
field has died down. Similar data is obtained for time step size 
corresponding to x = 2 for 50,000 time steps. As is evident 
from Fig. 7, the currents are stable for the entire duration. 

C. More challenging targets 

In this section, we present results for three challenging 
targets, a thin box, an even thinner NASA almond, and a cone 
sphere. All these scatterers are illuminated by a broadband 
pulse and analyzed for multiple transits across their surfaces. 
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Fig. 7: Current on two parallel plates 



1) Thin Box: Next, we analyze scattering from a thin box; 
again as mentioned earlier, this has been a challenge insofar 
as stability of the TDEFIE is concerned. To the best of our 
knowledge, stable results have been obtained only using exact 
integration and variations thereof [12, 13]. The dimension of 
the box are 0.5m x lm x 0.1m and is discretized with 
390 spatial unknowns. The box is excited by a field that is 
propagating along k = -y, polarized along u = z, and described 
by parameters f max = 211 MHz and f m { n - 1 MHz. The current 
is observed at (x,y,z) = (0.5,0,0) m, for 50,000 time steps 
for a time step size given by x - 0.5 and for 25,000 time 
steps for a time step size given by x — 1 • As is evident from 
Fig. 8, the two currents agree well with each other and they 
are stable. In addition to a full MOT solution, we conducted 
an eigen-analysis of the MOT system of equations, and as is 
evident from Fig. 9, all the eigenvalues lie within the unit 
circle. This again attests to the stability of the MOT system 
for this scatterer. 




Fig. 8: Current on thin box 

2) NASA Almond: Next, we analyze scattering from a thin 
almond; the almond fits in a box of dimension 2.17m x 
1.11m x 0.06m (aspect ratio 38:20:1), and is discretized using 
1140 unknowns. The target is illuminated by a field incident 
along k = z, polarized along u = -x, and characterized 
by fmax = 132 MHz and f m i n - 1 kHz. The current at 
(x,y,z) = (38.0,689,28.8) for 40,000 time steps with size 
determined by x — 1 is shown in Fig. 10. As is evident from 
this figure, the current does not show late time instability. 
The inset in Fig. 10 shows the features of the current until 
t = 03/is. Next, we extract RCS data at three different 
frequencies, in the x - z plane for £ [-180°, 180°], and 






Re( C) 



Fig. 9: Eigenvalues of thin box MOT matrix 



compare these results against similar data obtained using a 
frequency domain code. Again, as is evident in Fig. 11, the 
agreement between the two sets of data is excellent. 




15 



% s) 



Fig. 10: Current on NASA almond 
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Fig. 11: RCS of NASA almond 

3) Cone-Sphere: In this last example, we analyze scattering 
from a cone-sphere. The radius of the base of the cone is 
0.25m while its height is lm, and the scatterer is represented 
using 1008 spatial degrees of freedom. The incident field 
is propagating along k = z, is polarized along u = x, and 
is characterized by f max = 599 MHz and f m { n - 1 MHz. 
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The current observed at (x,y,z) = (0.329,0.335,0.125) m for 
40,000 time steps with the time step size corresponding to 
X = 1 is depicted in Fig. 12. As is evident from this figure, 
the currents exhibit late time stability for multiple transits 
across the geometry. An inset in Fig. 12 depicts features 
until t = 55ns. As before, RCS data in the x - z plane for 
6 £ [-180°, 180°] is extracted at three different frequencies 
and compared against similar data obtained using a frequency 
domain code. As is evident from Fig. 13, the agreement is 
excellent at all three frequencies. 




Fig. 12: Current on cone- sphere 




Fig. 13: RCS of cone-sphere 



V. Conclusion 

This paper presents a novel framework for constructing 
TDIEs; the crux of the approach presented in this paper lies 
in developing a separable spatio-temporal expansion for rep- 
resenting the convolution between the retarded potential and 
the source. As a result of this expansion, the discontinuities in 
the temporal basis set do not appear in the spatial integrands. 
This method, in concert with the correct variational form, has 
been applied to the analysis of scattering from number of 
challenging targets for multiple transits of the incident pulse 
across the object. The goal of this set of experiment was two- 
fold: (i) study stability behavior in early time, and (ii) analyze 
behavior when the incident field has completely died down. 
In all cases, the currents reach a DC floor, and remain there 
for the duration of the analysis; the DC floor is expected 
as it lies in the null space of the TDEFIE operator. These 
results demonstrate the viability of using this approach for 
TDIE analysis, and opens door to more challenging analysis. 



Extension of this approach to higher order geometries as well 
as integration with PWTD accelerators is underway and will 
be presented elsewhere. 

VI. Appendix 
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where [I] is the identity matrix. 
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